Chapter 4 Community composition
4.1 Taxonomy overview
4.1.1 Stacked barplot
# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
filter(count > 0) #filter 0 counts
# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)
ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ type, scales="free") + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
labs(fill="Phylum",y = "Relative abundance",x="Samples")+
scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label4.1.2 Phylum relative abundances
phylum_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
group_by(sample,phylum) %>%
summarise(relabun=sum(count))
phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| phylum | mean | sd |
|---|---|---|
| p__Bacteroidota | 0.376352935 | 0.205078492 |
| p__Bacillota_A | 0.250845368 | 0.165219604 |
| p__Bacillota | 0.116698680 | 0.146369286 |
| p__Pseudomonadota | 0.095767942 | 0.159976047 |
| p__Campylobacterota | 0.052429876 | 0.092626270 |
| p__Verrucomicrobiota | 0.029518627 | 0.072768229 |
| p__Desulfobacterota | 0.023304680 | 0.036356268 |
| p__Actinomycetota | 0.012917195 | 0.107880476 |
| p__Chlamydiota | 0.010449046 | 0.059349062 |
| p__Fusobacteriota | 0.010359534 | 0.028167393 |
| p__Cyanobacteriota | 0.009098787 | 0.016424861 |
| p__Bacillota_C | 0.004658039 | 0.006626014 |
| p__Spirochaetota | 0.003962784 | 0.012243051 |
| p__Bacillota_B | 0.002436629 | 0.004905803 |
| p__Elusimicrobiota | 0.001199879 | 0.006464719 |
phylum_arrange <- phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun)) %>%
arrange(-mean) %>%
select(phylum) %>%
pull()
phylum_summary %>%
filter(phylum %in% phylum_arrange) %>%
mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
geom_jitter(alpha=0.5) +
theme_minimal() +
theme(legend.position="none") +
labs(y="Phylum",x="Relative abundance")4.2 Taxonomy boxplot
4.2.1 Family
family_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,family) %>%
summarise(relabun=sum(count))
family_summary %>%
group_by(family) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| family | mean | sd |
|---|---|---|
| f__Bacteroidaceae | 2.189238e-01 | 0.1404298747 |
| f__Lachnospiraceae | 1.389745e-01 | 0.1089608538 |
| f__Tannerellaceae | 1.016053e-01 | 0.0802366634 |
| f__Helicobacteraceae | 5.197482e-02 | 0.0922005028 |
| f__Mycoplasmoidaceae | 3.651280e-02 | 0.0753015455 |
| f__Erysipelotrichaceae | 3.461428e-02 | 0.0450141949 |
| f__UBA3700 | 3.378612e-02 | 0.0555421417 |
| f__Marinifilaceae | 2.741603e-02 | 0.0272043836 |
| f__ | 2.740335e-02 | 0.0833966811 |
| f__DTU072 | 2.668734e-02 | 0.0526713763 |
| f__Rikenellaceae | 2.665357e-02 | 0.0462652505 |
| f__Enterobacteriaceae | 2.660196e-02 | 0.0913429878 |
| f__Coprobacillaceae | 2.521928e-02 | 0.0887342893 |
| f__Desulfovibrionaceae | 2.330468e-02 | 0.0363562684 |
| f__Ruminococcaceae | 1.838184e-02 | 0.0427420605 |
| f__LL51 | 1.738859e-02 | 0.0683714736 |
| f__Rhizobiaceae | 1.578130e-02 | 0.0766280394 |
| f__UBA3830 | 1.560351e-02 | 0.0451659123 |
| f__Mycobacteriaceae | 1.228390e-02 | 0.1079320546 |
| f__Akkermansiaceae | 1.213003e-02 | 0.0315247341 |
| f__Chlamydiaceae | 1.044905e-02 | 0.0593490616 |
| f__Fusobacteriaceae | 1.035953e-02 | 0.0281673926 |
| f__CAG-239 | 9.034415e-03 | 0.0150424841 |
| f__Enterococcaceae | 8.322037e-03 | 0.0463917038 |
| f__Gastranaerophilaceae | 7.638375e-03 | 0.0143794214 |
| f__Oscillospiraceae | 6.565519e-03 | 0.0077974910 |
| f__UBA1997 | 6.303806e-03 | 0.0307918484 |
| f__Streptococcaceae | 6.291980e-03 | 0.0340225424 |
| f__UBA1242 | 4.618700e-03 | 0.0152906152 |
| f__Brevinemataceae | 3.962784e-03 | 0.0122430505 |
| f__Acutalibacteraceae | 3.335081e-03 | 0.0108974502 |
| f__UBA660 | 3.159125e-03 | 0.0138778160 |
| f__Clostridiaceae | 2.808963e-03 | 0.0170355922 |
| f__RUG11792 | 2.784774e-03 | 0.0249664653 |
| f__Peptococcaceae | 2.436629e-03 | 0.0049058034 |
| f__MGBC116941 | 2.134904e-03 | 0.0093302936 |
| f__Acidaminococcaceae | 1.921070e-03 | 0.0050068934 |
| f__CAG-508 | 1.786832e-03 | 0.0063826540 |
| f__Anaerovoracaceae | 1.551174e-03 | 0.0036295746 |
| f__Moraxellaceae | 1.468042e-03 | 0.0096884882 |
| f__RUG14156 | 1.460412e-03 | 0.0045605184 |
| f__Staphylococcaceae | 1.345783e-03 | 0.0050584305 |
| f__Elusimicrobiaceae | 1.199879e-03 | 0.0064647185 |
| f__CAG-288 | 9.379861e-04 | 0.0059852367 |
| f__Anaerotignaceae | 8.884602e-04 | 0.0040242420 |
| f__CALVMC01 | 7.428930e-04 | 0.0043360033 |
| f__Eggerthellaceae | 6.332937e-04 | 0.0021152288 |
| f__Massilibacillaceae | 6.162148e-04 | 0.0016267915 |
| f__UBA1820 | 4.680638e-04 | 0.0012990360 |
| f__Arcobacteraceae | 4.550604e-04 | 0.0050029784 |
| f__CAG-274 | 4.466884e-04 | 0.0021903937 |
| f__Burkholderiaceae_C | 3.656162e-04 | 0.0047810524 |
| f__Muribaculaceae | 3.575580e-04 | 0.0009726557 |
| f__UBA932 | 3.379555e-04 | 0.0011520742 |
| f__Hepatoplasmataceae | 2.954146e-04 | 0.0038630477 |
| f__Rhodobacteraceae | 2.924483e-04 | 0.0038242573 |
| f__Weeksellaceae | 2.739210e-04 | 0.0031292133 |
| f__Eubacteriaceae | 1.627561e-04 | 0.0006691724 |
| f__Sphingobacteriaceae | 1.488164e-04 | 0.0012387571 |
| f__Devosiaceae | 1.472568e-04 | 0.0015006126 |
| f__Pumilibacteraceae | 1.262477e-04 | 0.0007602888 |
| f__WRAU01 | 9.491039e-05 | 0.0012411144 |
| f__Peptostreptococcaceae | 2.260586e-05 | 0.0002956099 |
family_arrange <- family_summary %>%
group_by(family) %>%
summarise(mean=sum(relabun)) %>%
arrange(-mean) %>%
select(family) %>%
pull()
family_summary %>%
left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
filter(family %in% family_arrange[1:20]) %>%
mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
scale_color_manual(values=phylum_colors[-8]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")4.2.2 Genus
genus_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,genus) %>%
summarise(relabun=sum(count)) %>%
filter(genus != "g__")
genus_summary %>%
group_by(genus) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| genus | mean | sd |
|---|---|---|
| g__Bacteroides | 1.337013e-01 | 0.0929327655 |
| g__Parabacteroides | 9.568074e-02 | 0.0803866834 |
| g__Phocaeicola | 6.854013e-02 | 0.0793558268 |
| g__Mycoplasmoides | 3.015602e-02 | 0.0749534641 |
| g__Helicobacter_J | 2.973811e-02 | 0.0592188466 |
| g__Odoribacter | 2.522456e-02 | 0.0267412147 |
| g__Roseburia | 2.356827e-02 | 0.0556733585 |
| g__NHYM01 | 2.223671e-02 | 0.0792620634 |
| g__Alistipes | 2.182523e-02 | 0.0282922451 |
| g__Coprobacillus | 1.990305e-02 | 0.0873936500 |
| g__Agrobacterium | 1.578130e-02 | 0.0766280394 |
| g__Corynebacterium | 1.228390e-02 | 0.1079320546 |
| g__Akkermansia | 1.213003e-02 | 0.0315247341 |
| g__Fusobacterium_A | 1.026752e-02 | 0.0281722895 |
| g__Kineothrix | 8.698172e-03 | 0.0406735379 |
| g__Proteus | 8.556721e-03 | 0.0677872967 |
| g__Dielma | 8.477779e-03 | 0.0090773537 |
| g__CAG-95 | 8.014504e-03 | 0.0203997532 |
| g__JAAYNV01 | 7.901098e-03 | 0.0194443338 |
| g__UBA866 | 7.175409e-03 | 0.0291873132 |
| g__Desulfovibrio | 6.936071e-03 | 0.0210368286 |
| g__Enterococcus | 6.919414e-03 | 0.0453969341 |
| g__Ureaplasma | 6.356785e-03 | 0.0137373135 |
| g__Lactococcus | 6.291980e-03 | 0.0340225424 |
| g__Lacrimispora | 5.993828e-03 | 0.0097829951 |
| g__Parabacteroides_B | 5.924566e-03 | 0.0099792665 |
| g__CALXRO01 | 5.698293e-03 | 0.0306767185 |
| g__Citrobacter | 5.650169e-03 | 0.0332631355 |
| g__NSJ-61 | 5.495080e-03 | 0.0197757124 |
| g__Breznakia | 5.440147e-03 | 0.0235692957 |
| g__Clostridium_AQ | 5.263895e-03 | 0.0121113284 |
| g__Salmonella | 5.073541e-03 | 0.0183551010 |
| g__Bilophila | 4.925550e-03 | 0.0088098210 |
| g__Hungatella_A | 4.755524e-03 | 0.0095127458 |
| g__MGBC136627 | 4.313746e-03 | 0.0162110012 |
| g__Escherichia | 4.139378e-03 | 0.0264569200 |
| g__UMGS1251 | 4.111189e-03 | 0.0072426766 |
| g__Hungatella | 4.068243e-03 | 0.0189714607 |
| g__Clostridium_Q | 3.965072e-03 | 0.0052000165 |
| g__Brevinema | 3.962784e-03 | 0.0122430505 |
| g__Thomasclavelia | 3.856935e-03 | 0.0108480346 |
| g__Scatousia | 3.607252e-03 | 0.0102197530 |
| g__Enterocloster | 3.578040e-03 | 0.0047221634 |
| g__Mailhella | 3.569833e-03 | 0.0101940649 |
| g__Copromonas | 3.557270e-03 | 0.0050035176 |
| g__Ventrimonas | 3.472264e-03 | 0.0070914148 |
| g__Caccovivens | 3.304213e-03 | 0.0121527504 |
| g__Lawsonia | 3.250547e-03 | 0.0116543597 |
| g__Fournierella | 3.179598e-03 | 0.0062039376 |
| g__Limenecus | 3.126281e-03 | 0.0065432982 |
| g__MGBC133411 | 3.006834e-03 | 0.0074209199 |
| g__Mucinivorans | 2.866176e-03 | 0.0371005376 |
| g__Sarcina | 2.808963e-03 | 0.0170355922 |
| g__Acetatifactor | 2.706113e-03 | 0.0058179611 |
| g__Eisenbergiella | 2.665864e-03 | 0.0067990748 |
| g__Bacteroides_G | 2.651345e-03 | 0.0344120334 |
| g__CAJLXD01 | 2.603014e-03 | 0.0087025932 |
| g__Blautia | 2.528781e-03 | 0.0061107674 |
| g__C-19 | 2.248901e-03 | 0.0048466847 |
| g__Butyricimonas | 2.191467e-03 | 0.0049842986 |
| g__Velocimicrobium | 2.175479e-03 | 0.0066412559 |
| g__CAZU01 | 2.170696e-03 | 0.0065598332 |
| g__MGBC116941 | 2.134904e-03 | 0.0093302936 |
| g__Intestinimonas | 2.057548e-03 | 0.0039214573 |
| g__Negativibacillus | 2.044877e-03 | 0.0054857591 |
| g__Rikenella | 1.962168e-03 | 0.0036910712 |
| g__Phascolarctobacterium | 1.921070e-03 | 0.0050068934 |
| g__RGIG6463 | 1.768909e-03 | 0.0039495868 |
| g__JALFVM01 | 1.721982e-03 | 0.0038441918 |
| g__Oscillibacter | 1.492792e-03 | 0.0024898538 |
| g__Pseudoflavonifractor | 1.485397e-03 | 0.0027590442 |
| g__Acinetobacter | 1.468042e-03 | 0.0096884882 |
| g__Citrobacter_A | 1.376632e-03 | 0.0060011629 |
| g__Staphylococcus | 1.345783e-03 | 0.0050584305 |
| g__RGIG4733 | 1.288541e-03 | 0.0040266511 |
| g__UBA1436 | 1.199879e-03 | 0.0064647185 |
| g__Lachnotalea | 1.194066e-03 | 0.0048896855 |
| g__Ruthenibacterium | 1.187962e-03 | 0.0053331515 |
| g__14-2 | 1.170787e-03 | 0.0095739735 |
| g__Beduini | 1.160220e-03 | 0.0025025910 |
| g__Scatocola | 1.108079e-03 | 0.0044726561 |
| g__Faecisoma | 1.072888e-03 | 0.0055281036 |
| g__Enterococcus_A | 1.071313e-03 | 0.0098572488 |
| g__Faecimonas | 9.764964e-04 | 0.0082595968 |
| g__RGIG9287 | 9.526092e-04 | 0.0092684203 |
| g__CAG-345 | 9.379861e-04 | 0.0059852367 |
| g__Blautia_A | 9.100332e-04 | 0.0028927770 |
| g__RGIG1896 | 8.246024e-04 | 0.0062408705 |
| g__CAG-269 | 7.889086e-04 | 0.0046906408 |
| g__Marseille-P3106 | 7.845423e-04 | 0.0017535790 |
| g__WRHT01 | 6.354364e-04 | 0.0026829604 |
| g__Eggerthella | 6.332937e-04 | 0.0021152288 |
| g__IOR16 | 6.324101e-04 | 0.0020540977 |
| g__Anaerotruncus | 6.192165e-04 | 0.0016862064 |
| g__RUG14156 | 6.079829e-04 | 0.0022026464 |
| g__CHH4-2 | 6.073170e-04 | 0.0019890673 |
| g__Serratia_A | 5.792071e-04 | 0.0075741155 |
| g__CAG-56 | 4.857941e-04 | 0.0016254203 |
| g__Merdimorpha | 4.680638e-04 | 0.0012990360 |
| g__MGBC140009 | 4.624605e-04 | 0.0023930648 |
| g__CALURL01 | 4.581073e-04 | 0.0016646244 |
| g__Aliarcobacter | 4.550604e-04 | 0.0050029784 |
| g__Scatenecus | 4.467347e-04 | 0.0019650337 |
| g__RGIG8482 | 4.347613e-04 | 0.0029582243 |
| g__Enterobacter | 4.025794e-04 | 0.0041076314 |
| g__Klebsiella | 4.007019e-04 | 0.0048624260 |
| g__Caccenecus | 3.895102e-04 | 0.0017702443 |
| g__Alcaligenes | 3.656162e-04 | 0.0047810524 |
| g__Plesiomonas | 3.590755e-04 | 0.0026947991 |
| g__HGM05232 | 3.575580e-04 | 0.0009726557 |
| g__JAHHSE01 | 3.550107e-04 | 0.0014818046 |
| g__Egerieousia | 3.379555e-04 | 0.0011520742 |
| g__Emergencia | 3.373704e-04 | 0.0017351293 |
| g__Enterococcus_B | 3.313107e-04 | 0.0022138492 |
| g__Stoquefichus | 2.990680e-04 | 0.0020385612 |
| g__Hepatoplasma | 2.954146e-04 | 0.0038630477 |
| g__Paracoccus | 2.924483e-04 | 0.0038242573 |
| g__Moheibacter | 2.739210e-04 | 0.0031292133 |
| g__Scatomorpha | 2.610126e-04 | 0.0010128258 |
| g__UBA7185 | 2.405856e-04 | 0.0014474682 |
| g__Eubacterium | 1.627561e-04 | 0.0006691724 |
| g__Sphingobacterium | 1.488164e-04 | 0.0012387571 |
| g__Devosia | 1.472568e-04 | 0.0015006126 |
| g__Anaerosporobacter | 1.437105e-04 | 0.0012673026 |
| g__Caccomorpha | 1.366945e-04 | 0.0010479478 |
| g__UBA2658 | 1.292159e-04 | 0.0007163845 |
| g__Protoclostridium | 1.262477e-04 | 0.0007602888 |
| g__Angelakisella | 1.253643e-04 | 0.0009168117 |
| g__Cetobacterium_A | 9.201325e-05 | 0.0008667719 |
| g__Rahnella | 6.395024e-05 | 0.0008362579 |
| g__Peptostreptococcus | 2.260586e-05 | 0.0002956099 |
genus_arrange <- genus_summary %>%
group_by(genus) %>%
summarise(mean=sum(relabun)) %>%
filter(genus != "g__")%>%
arrange(-mean) %>%
select(genus) %>%
mutate(genus= sub("^g__", "", genus)) %>%
pull()
genus_summary %>%
left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
mutate(genus= sub("^g__", "", genus)) %>%
filter(genus %in% genus_arrange[1:20]) %>%
mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")4.3 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))4.3.1 Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric ) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.3.2 Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.3.3 Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.3.4 Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.3.5 Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.3.6 Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.3.7 Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)4.4 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)4.5 Permanovas
4.5.1 1. Are the wild populations similar?
4.5.1.1 Wild: P.muralis vs P.liolepis
wild <- meta %>%
filter(time_point == "0_Wild")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))
wild_nmds <- sample_metadata %>%
filter(time_point == "0_Wild")4.5.1.3 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.631953 | 0.207198 | 6.795072 | 0.001 |
| Residual | 26 | 6.244347 | 0.792802 | NA | NA |
| Total | 27 | 7.876300 | 1.000000 | NA | NA |
4.5.1.4 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 2.018797 | 0.2566342 | 8.976049 | 0.001 |
| Residual | 26 | 5.847641 | 0.7433658 | NA | NA |
| Total | 27 | 7.866438 | 1.0000000 | NA | NA |
4.5.2 2. Do the antibiotics work?
4.5.2.1 Acclimation vs antibiotics
treat <- meta %>% #antibiotics samples giving error when plotting
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))
treat_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")4.5.2.2 Number of samples used
[1] 55
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts_func, q=1, dist=dist)4.5.2.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.042510 | 0.0911466 | 6.563244 | 0.001 |
| species | 1 | 2.300885 | 0.1026765 | 7.393486 | 0.001 |
| individual | 26 | 9.974365 | 0.4451038 | 1.232725 | 0.003 |
| Residual | 26 | 8.091314 | 0.3610731 | NA | NA |
| Total | 54 | 22.409075 | 1.0000000 | NA | NA |
4.5.2.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.189257 | 0.1017464 | 8.979451 | 0.001 |
| species | 1 | 3.046568 | 0.1415902 | 12.495800 | 0.001 |
| individual | 26 | 9.941985 | 0.4620568 | 1.568386 | 0.001 |
| Residual | 26 | 6.338992 | 0.2946066 | NA | NA |
| Total | 54 | 21.516802 | 1.0000000 | NA | NA |
4.5.2.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.1642230 | 0.2000313 | 18.460841 | 0.001 |
| species | 1 | 0.8686306 | 0.0802844 | 7.409427 | 0.001 |
| individual | 26 | 4.7385071 | 0.4379630 | 1.554596 | 0.001 |
| Residual | 26 | 3.0480626 | 0.2817214 | NA | NA |
| Total | 54 | 10.8194233 | 1.0000000 | NA | NA |
4.5.2.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.2832796 | 0.3550091 | 33.1502018 | 0.001 |
| species | 1 | 0.0227598 | 0.0035387 | 0.3304425 | 0.708 |
| individual | 26 | 2.3347730 | 0.3630154 | 1.3037622 | 0.191 |
| Residual | 26 | 1.7907967 | 0.2784368 | NA | NA |
| Total | 54 | 6.4316092 | 1.0000000 | NA | NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))4.5.3 3. Do the FMT work?
4.5.3.1 Comparison between FMT2 vs Post-FMT2
transplant3 <- meta %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")
transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))
transplant3.counts_func <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts_func)),sort(as.character(rownames(transplant3))))
transplant3_nmds <- sample_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")4.5.3.2 Number of samples used
[1] 54
beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts_func, q=1, dist=dist)4.5.3.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9748558 | 0.0733017 | 5.103443 | 0.001 |
| time_point | 1 | 0.7043146 | 0.0529590 | 3.687140 | 0.001 |
| type | 1 | 1.4084037 | 0.1059011 | 7.373099 | 0.001 |
| individual | 25 | 7.1553504 | 0.5380273 | 1.498352 | 0.001 |
| Residual | 16 | 3.0563076 | 0.2298108 | NA | NA |
| Total | 44 | 13.2992322 | 1.0000000 | NA | NA |
4.5.3.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.1035252 | 0.0842380 | 6.734701 | 0.001 |
| time_point | 1 | 0.7327432 | 0.0559342 | 4.471856 | 0.001 |
| type | 1 | 1.8565350 | 0.1417193 | 11.330241 | 0.001 |
| individual | 25 | 6.7855728 | 0.5179794 | 1.656466 | 0.001 |
| Residual | 16 | 2.6217058 | 0.2001290 | NA | NA |
| Total | 44 | 13.1000821 | 1.0000000 | NA | NA |
4.5.3.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1387836 | 0.0608025 | 2.904615 | 0.126 |
| time_point | 1 | 0.1416869 | 0.0620745 | 2.965380 | 0.002 |
| type | 1 | 0.1763572 | 0.0772639 | 3.690999 | 0.002 |
| individual | 25 | 1.0612162 | 0.4649298 | 0.888412 | 0.393 |
| Residual | 16 | 0.7644858 | 0.3349292 | NA | NA |
| Total | 44 | 2.2825297 | 1.0000000 | NA | NA |
4.5.3.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0003437 | -0.0001753 | -0.0062959 | 0.903 |
| time_point | 1 | 0.0253408 | 0.0129285 | 0.4642107 | 0.713 |
| type | 1 | 0.1178144 | 0.0601071 | 2.1582044 | 0.137 |
| individual | 25 | 0.9438371 | 0.4815314 | 0.6915940 | 0.747 |
| Residual | 16 | 0.8734253 | 0.4456084 | NA | NA |
| Total | 44 | 1.9600740 | 1.0000000 | NA | NA |
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylo_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_func_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")4.5.4 4. Are there differences between the control and the treatment group?
4.5.4.1 Post-FMT2
post2 <- meta %>%
filter(time_point == "6_Post-FMT2")
post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))
post2.counts_func <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts_func)),sort(as.character(rownames(post2))))
post2_nmds <- sample_metadata %>%
filter(time_point == "6_Post-FMT2")4.5.4.2 Number of samples used
[1] 28
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts_func, q=1, dist=dist)4.5.4.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.8966107 | 0.1132233 | 3.515588 | 0.001 |
| type | 1 | 0.6463814 | 0.0816245 | 2.534445 | 0.001 |
| Residual | 25 | 6.3759661 | 0.8051521 | NA | NA |
| Total | 27 | 7.9189582 | 1.0000000 | NA | NA |
4.5.4.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9398968 | 0.1224370 | 4.112304 | 0.001 |
| type | 1 | 1.0227481 | 0.1332297 | 4.474800 | 0.001 |
| Residual | 25 | 5.7139316 | 0.7443333 | NA | NA |
| Total | 27 | 7.6765766 | 1.0000000 | NA | NA |
4.5.4.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1200820 | 0.1454923 | 4.647187 | 0.001 |
| type | 1 | 0.0592745 | 0.0718175 | 2.293931 | 0.031 |
| Residual | 25 | 0.6459931 | 0.7826902 | NA | NA |
| Total | 27 | 0.8253496 | 1.0000000 | NA | NA |
4.5.4.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0116738 | 0.0168395 | 0.4229082 | 0.530 |
| type | 1 | -0.0085273 | -0.0123007 | -0.3089208 | 0.915 |
| Residual | 25 | 0.6900904 | 0.9954612 | NA | NA |
| Total | 27 | 0.6932368 | 1.0000000 | NA | NA |
p0<-beta_richness_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")